//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry             : env klmperry
global storageb				: env storageb

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"
global dataihdp         = "$storageb/dc_data/"

cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// low birthweight
rename tot_siblings_natural siblings_natural

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

rename cum_avg_daycare_36m_sum c
rename parenting_ages13_std    p

global inputs_std c p

tab   site, gen(site_)
global baseline_child      sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare siblings_natural employed_adult 

// discretize
drop    bw 
gen     bw = 0 if bwg == "L"
replace bw = 1 if bwg == "H"

foreach var of varlist bw anga mage meduc  {
	summ    `var', d
	gen     `var'_p = 0 if `var' <= r(p50) 
	replace `var'_p = 1 if `var' >  r(p50) & `var' !=.
	drop    `var'
	rename  `var'_p `var'
}
replace employed_adult   = 1 if employed_adult   != 0
replace siblings_natural = 1 if siblings_natural >  1

keep  tg twin $baseline_child $baseline_mother $baseline_household $outputs  site_1-site_8

foreach var of varlist twin $baseline_child $baseline_mother $baseline_household  site_1-site_8 {
	gen `var'_tg = `var'*tg
}

// factors
foreach var of varlist $baseline_child $baseline_mother $baseline_household {
	summ  `var' 
	gen   `var'_std = (`var' - r(mean))/r(sd)
}
sem (_cons@0 X -> sex_std) (_cons@0 X -> bw_std) (_cons@0 X -> anga_std) (_cons@0 X -> black_std) (_cons@0 X -> hispanic_std)
predict factor_child, latent
sem (_cons@0 X -> mage_std) (_cons@0 X -> meduc_std) (_cons@0 X -> works_std) (_cons@0 X -> married_std)
predict factor_mother, latent
sem (_cons@0 X -> welfare_std) (_cons@0 X -> siblings_natural_std) (_cons@0 X -> employed_adult_std)
predict factor_household, latent
foreach var in child mother household {
	gen factor_`var'_tg = factor_`var'*tg
}

// also average
egen iq3 = rowmean(iqcage stndscor)

// by child
foreach vary of varlist iqcage stndscor iq3 {
	matrix child_`vary' = J(1,6,.)
	foreach var of varlist twin $baseline_child {
		# delimit
		bootstrap, reps($bootstraps): reg `vary' tg twin $baseline_child factor_mother factor_household 
									  factor_mother_tg factor_household_tg `var'_tg;
		# delimit cr
		
		// b, V
		matrix b = e(b)
		matrix V = e(V)
		
		// tg
		matrix se`var'  = sqrt(V[1,1])
		matrix md`var'  = b[1,1]
		matrix mdt`var' = abs(md`var'[1,1]/se`var'[1,1])
		matrix mdp`var' = 2*(1 - normal(mdt`var'[1,1]))
		
		// tg * var
		matrix sei`var'  = sqrt(V[12,12])
		matrix mi`var'  = b[1,12]
		matrix mit`var' = abs(mi`var'[1,1]/sei`var'[1,1])
		matrix mip`var' = 2*(1 - normal(mit`var'[1,1]))
		
		// ratio
		nlcom _b[`var'_tg]/_b[tg]
		matrix ser`var' = sqrt(r(V)[1,1])
		matrix mr`var'  = abs(r(b)[1,1])
		matrix mrt`var' = abs(mr`var'[1,1]/ser`var'[1,1])
		matrix mrp`var' = 2*(1 - normal(mrt`var'[1,1])) 
		
		matrix child_`vary' = [child_`vary' \ [md`var',mdp`var',mi`var',mip`var',mr`var',mrp`var']]
	}
	matrix colnames child_`vary' = mdc_`vary' mdpc_`vary' mic_`vary' mipc_`vary' mrc_`vary' mrpc_`vary'

	// by mother
	matrix mother_`vary' = J(1,6,.)
	foreach var of varlist twin $baseline_mother {
		# delimit
		bootstrap, reps($bootstraps): reg `vary' tg twin $baseline_mother factor_child factor_household 
									  factor_child_tg factor_household_tg `var'_tg;
		# delimit cr
		
		// b, V
		matrix b = e(b)
		matrix V = e(V)
		
		// tg
		matrix se`var'  = sqrt(V[1,1])
		matrix md`var'  = b[1,1]
		matrix mdt`var' = abs(md`var'[1,1]/se`var'[1,1])
		matrix mdp`var' = 2*(1 - normal(mdt`var'[1,1]))
		
		// tg * var
		matrix sei`var'  = sqrt(V[11,11])
		matrix mi`var'  = b[1,11]
		matrix mit`var' = abs(mi`var'[1,1]/sei`var'[1,1])
		matrix mip`var' = 2*(1 - normal(mit`var'[1,1]))
		
		// ratio
		nlcom _b[`var'_tg]/_b[tg]
		matrix ser`var' = sqrt(r(V)[1,1])
		matrix mr`var'  = abs(r(b)[1,1])
		matrix mrt`var' = abs(mr`var'[1,1]/ser`var'[1,1])
		matrix mrp`var' = 2*(1 - normal(mrt`var'[1,1])) 
		
		matrix mother_`vary' = [mother_`vary' \ [md`var',mdp`var',mi`var',mip`var',mr`var',mrp`var']]
	}
	matrix colnames mother_`vary' = mdm_`vary' mdpm_`vary' mim_`vary' mipm_`vary' mrm_`vary' mrpm_`vary'

	// by household
	matrix household_`vary' = J(1,6,.)
	foreach var of varlist twin $baseline_household {
		# delimit
		bootstrap, reps($bootstraps): reg `vary' tg twin $baseline_household factor_child factor_mother 
									  factor_child_tg factor_mother_tg `var'_tg;
		# delimit cr
		
		// b, V
		matrix b = e(b)
		matrix V = e(V)
		
		// tg
		matrix se`var'  = sqrt(V[1,1])
		matrix md`var'  = b[1,1]
		matrix mdt`var' = abs(md`var'[1,1]/se`var'[1,1])
		matrix mdp`var' = 2*(1 - normal(mdt`var'[1,1]))
		
		// tg * var
		matrix sei`var'  = sqrt(V[10,10])
		matrix mi`var'  = b[1,10]
		matrix mit`var' = abs(mi`var'[1,1]/sei`var'[1,1])
		matrix mip`var' = 2*(1 - normal(mit`var'[1,1]))
		
			
		// ratio
		nlcom _b[`var'_tg]/_b[tg]
		matrix ser`var' = sqrt(r(V)[1,1])
		matrix mr`var'  = abs(r(b)[1,1])
		matrix mrt`var' = abs(mr`var'[1,1]/ser`var'[1,1])
		matrix mrp`var' = 2*(1 - normal(mrt`var'[1,1])) 
		
		matrix household_`vary' = [household_`vary' \ [md`var',mdp`var',mi`var',mip`var',mr`var',mrp`var']]
	}
	matrix colnames household_`vary' = mdh_`vary' mdph_`vary' mih_`vary' miph_`vary' mrh_`vary' mrph_`vary'
}

// graph
clear
foreach vary in iqcage stndscor iq3 {
	foreach var in child mother household {
		matrix `var'_`vary' = `var'_`vary'[2...,1...]
		svmat `var'_`vary', names(col)
	}
}
gen     nd = _n*3 - 2
gen     ni = nd + 1
egen    ndni = rowmean(nd ni)

// generate ratio
foreach vary in iqcage stndscor iq3 {
	foreach var in c m h { 
		tostring  mr`var'_`vary', gen(mr`var'str_`vary') force format(%15.2fc) 
		replace   mr`var'str_`vary' = "|{&Delta}|/|{&Delta}{sup:Int}| = " + mr`var'str_`vary'
		
		tostring  mrp`var'_`vary', gen(mrp`var'str_`vary') force format(%15.2fc) 
		replace   mrp`var'str_`vary' = "{it:p} = " + mrp`var'str_`vary'
	}
}
gen cons1 = -11
gen cons2 = -12

// child
cd $output
foreach vary in iqcage stndscor iq3 {
	#delimit
		   twoway (bar     mdc_`vary' nd if _n <= 6,  color(gs10) 				barwidth(.9)                                         yline(0  , lpattern(solid) lwidth(vthick) lcolor(gs14)))	   
				  (scatter mdc_`vary' nd if mdpc_`vary'  < 0.05                   & _n <= 6, msize(medlarge)  mcolor(gs0) msymbol(circle)  yline(12 , lpattern(solid) lcolor(gs14)))
				  (scatter mdc_`vary' nd if mdpc_`vary'  >= .05 & mdpc_`vary'    <= 0.10 & _n <= 6, msize(medlarge)  mcolor(gs0) msymbol(square)  yline(-12, lpattern(solid) lcolor(gs14)))
				  
				  (bar     mic_`vary' ni if _n <= 6,  fcolor(none) lcolor(gs0) 				barwidth(.9))	   
				  (scatter mic_`vary' ni if mipc_`vary'  < 0.05                   & _n <= 6, msize(medlarge)  mcolor(gs0) msymbol(circle))
				  (scatter mic_`vary' ni if mipc_`vary'  >= .05 & mipc_`vary'    <= 0.10 & _n <= 6, msize(medlarge)  mcolor(gs0) msymbol(square))
				  
				  (scatter cons1 ndni    if _n <= 6, msymbol(none) mlabel(mrcstr_`vary')  mlabposition(12) mlabcolor(gs0) mlabsize(small))
				  (scatter cons2 ndni    if _n <= 6, msymbol(none) mlabel(mrpcstr_`vary') mlabposition(12) mlabcolor(gs0) mlabsize(small))
				  
				  ,
		   
			  legend(label(1 "{&Delta}: Overall Treatment Effect") label(2 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
					 label(3 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
					 label(4 "{&Delta}{sup:Int}: Interaction Treatment Effect")
					 cols(2) rows(2) order(1 4 2 3) size(small) position(6) region(lcolor(gs0)))
					 
			  ylabel(-12[6]12, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
			  xlabel(1.5 "Twin" 4.5 "Male" 7.5 `" "Birth" "Weight" "' 10.5 `" "Gestational" "Age" "' 13.5 "Black" 16.5 "Hispanic", labsize(small) grid glcolor(gs14) glpattern(solid))
			  ytitle("Treatment Effect", size(medlarge))
			  xtitle("")
			  graphregion(color(white)) plotregion(fcolor(white)); 
	#delimit cr
	graph export child_mlmot_`vary'.eps, replace
}

// mother
foreach vary in iqcage stndscor iq3 {
	#delimit
		   twoway (bar     mdm_`vary' nd if _n <= 5,  color(gs10) 				barwidth(.9)                                         yline(0  , lpattern(solid) lwidth(vthick) lcolor(gs14)))	   
				  (scatter mdm_`vary' nd if mdpc_`vary'  < 0.05                   & _n <= 5, msize(medlarge)  mcolor(gs0) msymbol(circle)  yline(12 , lpattern(solid) lcolor(gs14)))
				  (scatter mdm_`vary' nd if mdpc_`vary'  >= .05 & mdpm_`vary'    <= 0.10 & _n <= 5, msize(medlarge)  mcolor(gs0) msymbol(square)  yline(-12, lpattern(solid) lcolor(gs14)))
				  
				  (bar     mim_`vary' ni if _n <= 5,  fcolor(none) lcolor(gs0) 				barwidth(.9))	   
				  (scatter mim_`vary' ni if mipc_`vary'  < 0.05                   & _n <= 5, msize(medlarge)  mcolor(gs0) msymbol(circle))
				  (scatter mim_`vary' ni if mipc_`vary'  >= .05 & mipm_`vary'    <= 0.10 & _n <= 5, msize(medlarge)  mcolor(gs0) msymbol(square))
				  
				  (scatter cons1 ndni    if _n <= 5, msymbol(none) mlabel(mrmstr_`vary')  mlabposition(12) mlabcolor(gs0) mlabsize(small))
				  (scatter cons2 ndni    if _n <= 5, msymbol(none) mlabel(mrpmstr_`vary') mlabposition(12) mlabcolor(gs0) mlabsize(small))
				  
				  ,
		   
			  legend(label(1 "{&Delta}: Overall Treatment Effect") label(2 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
					 label(3 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
					 label(4 "{&Delta}{sup:Int}: Interaction Treatment Effect")
					 cols(2) rows(2) order(1 4 2 3) size(small) position(6) region(lcolor(gs0)))
					 
			  ylabel(-12[6]12, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
			  xlabel(1.5 "Twin" 4.5 "Age" 7.5 "Education" 10.5 "Works" 13.5 "Married", labsize(small) grid glcolor(gs14) glpattern(solid))
			  ytitle("Treatment Effect", size(medlarge))
			  xtitle("")
			  graphregion(color(white)) plotregion(fcolor(white)); 
	#delimit cr
	graph export mother_mlmot_`vary'.eps, replace
}

// household
foreach vary in iqcage stndscor iq3 {
	#delimit
		   twoway (bar     mdh_`vary' nd if _n <= 4,  color(gs10) 				barwidth(.9)                                         yline(0  , lpattern(solid) lwidth(vthick) lcolor(gs14)))	   
				  (scatter mdh_`vary' nd if mdph_`vary'  < 0.05                   & _n <= 4, msize(medlarge)  mcolor(gs0) msymbol(circle)  yline(12 , lpattern(solid) lcolor(gs14)))
				  (scatter mdh_`vary' nd if mdph_`vary'  >= .05 & mdph_`vary'    <= 0.10 & _n <= 4, msize(medlarge)  mcolor(gs0) msymbol(square)  yline(-12, lpattern(solid) lcolor(gs14)))
				  
				  (bar     mih_`vary' ni if _n <= 4,  fcolor(none) lcolor(gs0) 				barwidth(.9))	   
				  (scatter mih_`vary' ni if miph_`vary'  < 0.05                   & _n <= 4, msize(medlarge)  mcolor(gs0) msymbol(circle))
				  (scatter mih_`vary' ni if miph_`vary'  >= .05 & miph_`vary'    <= 0.10 & _n <= 4, msize(medlarge)  mcolor(gs0) msymbol(square))
				  
				  (scatter cons1 ndni    if _n <= 4, msymbol(none) mlabel(mrhstr_`vary')  mlabposition(12) mlabcolor(gs0) mlabsize(small))
				  (scatter cons2 ndni    if _n <= 4, msymbol(none) mlabel(mrphstr_`vary') mlabposition(12) mlabcolor(gs0) mlabsize(small))
				  
				  ,
		   
			  legend(label(1 "{&Delta}: Overall Treatment Effect") label(2 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
					 label(3 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
					 label(4 "{&Delta}{sup:Int}: Interaction Treatment Effect")
					 cols(2) rows(2) order(1 4 2 3) size(small) position(6) region(lcolor(gs0)))
					 
			  ylabel(-12[6]12, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
			  xlabel(1.5 "Twin" 4.5  `" "Social" "Programs" "' 7.5 "Siblings" 10.5  `" "Employed" "Adults" "' 
								, labsize(small) grid glcolor(gs14) glpattern(solid))
			  ytitle("Treatment Effect", size(medlarge))
			  xtitle("")
			  graphregion(color(white)) plotregion(fcolor(white)); 
	#delimit cr
	graph export household_mlmot_`vary'.eps, replace
}
